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Abstract 


One of the areas of great emphasis in Bhabha Atomic Research Center (BARC) 
Mumbai is dropwise condensation of metal vapor on the lower side of the cooled 
substrate. In the present work, drop shapes are generated on and below the 
inclined plate. This is done for 4 different angles of inclination and 3 different 
liquids of different volumes. In the static case, stability analysis for the drop 
above and below the inclined plate is performed. Advancing and receding angles 
are calculated for the drop on and below the inclined plate. Center of gravity 
and moment calculations are performed for both configurations. 

For the dynamic analysis of the drop movement, a semi-circular drop which is 
kept on a fiat plate moving with uniform velocity is first modeled. Pressure and 
velocity values are calculated for the semi-circular drop, by solving the governing 
equations with appropriate boundary conditions. Grid independence test and 
code validation has also been carried out for the semi-circular drop. Coefficient 
of friction between the drop and plate is calculated. The assumptions taken for 
these calculations are that the flow is steady, incompressible and the Reynolds 
number is small enough for the creeping flow approximation to be valid. 

The next step is grid generation for handling deformed drops. This approach 
is required when we have a complex geometry in the physical domain. We can 
transform the complex geometry from the physical to the computational domain 
by mapping each grid point in the physical domain to a convenient grid point in 
the computational domain. After the solution is obtained for all grid points in the 
computational domain, the solution is again mapped into the physical domain 
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through the transformation relations. The governing equations and boundary 
conditions are also transformed into the computational domain. In the present 
work, the physical domain is a deformed drop geometry while the computational 
domain is a semi-circular drop. Coefficient of friction has been calculated on 
the plate for the deformed drop in the physical domain. Results show that the 
coefficient of friction increases as Reynolds number and Weber number decreases. 
The effect of deformation is seen in the local pressure and i^elocity fields, but the 
overall impact on the coefficient of friction is small. 

A proposal for solving high Reynolds number flows, energy equation, treat- 
ment of a cluster of drops and coalescence of adjacent drops has also been pre- 
sented. 
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r* Radial distance, m 
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Chapter 1 
Introduction 


Condensation of metallic vapor on the lower side of a cooled substrate, is one 
of the important areas of research in Bhabha Atomic Research Center (BARC), 
Mumbai. The process of vaporization and condensation takes place in a vacuum 
chamber, this is because the electron beam used to vaporize the metal from the 
crucible, does not interacts with any other material present, and thus efficiency 
achieved in vaporization is very high. In condensation, the condensation sites are 
either available on the substrate or on the dust particles suspended in the air, but 
as there is vacuum the process of condensation takes principally on the substrate. 



Figure 1.1: Typical arrangement involving condensation of metal vapor. 



2 


Arrangement involving condensation of metal vapor is shown in Figure 1.1. 
The metal material to be vaporized is kept in a water cooled crucible. The crucible 
material should be such that it should have limited or no solubility wdth the metal 
kept in it. The electron beam used to vaporize the metal from the crucible has a 
power rating of about 3kW and a peak value of 30kW. The electron beam, which 
is incident on the metal material in the crucible, makes an angle of 45° to the 
horizontal. Crucible is used because it facilitates the batch evaporation of the 
metal for a long period of time. 

The stagnation temperature of the vapor is 1500K and the temperature of 
the substrate is 900K. So there is a large temperature gradient due to which 
condensation on the substrate is independent of the temperature difference. As 
the metal vapors rise in the vacuum chamber and reach the substrate, they are 
condensed on it. 

The substrate should have such surface composition and thermo-physical 
properties that it permits dropwise condensation of metal vapors on it. It is tilted 
to 7 — 15° to the horizontal to make the drainage of drops effective. Although the 
distribution of condensation sites is usually random, on an average, the density 
of nucleation sites is of the order of lO^sites/cm^. However the substrate surface 
can be made to have a preferred orientation. 

When a drop is condensed on the substrate, its initial growth is due to the 
additional condensation of the vapor. But as time passes the drop size increases 
and inter-drop distance decreases. At a particular value of the drop radius and 
inter-drop distance the drop wdll coalesce with one or more drops around it. Con- 
densation results in sudden increase in the radius of the drop and also its surface 
area. As the sum of the surface areas occupied by the drops individually before 
coalescence is greater than the surface area occupied by the single drop after co- 
alescence, so the condensation is accelerated as now more surface is available for 
condensation. 

After the process of condensation and coalescence have been initiated, at any 
instant of time, the substrate contains drops of various sizes. When a critical 
size of the drops is achieved, it can either fall-off or come down the substrate. 
When a drop comes down the substrate, it can slide or roll-down. But sliding 
motion is more prevalent than rolling motion in practical applications. As the 
drops slides down the substrate, they are collected in a receptor. Moreover, the 
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sliding process can further facilitate the coalescence. 

Dropwise condensation is a complex process. In order to simulate the whole 
process from dropwise condensation till the drops are collected in a receptor, it 
is necessary to develop a model for an individual drop. Thus, the present work 
deals with the modeling of an individual drop. 

This work starts from the simulation of drop shapes on an inclined surface, 
both above and below^ Static drops involve estimation of threshold angle of 
inclination of the flat plate for which the drop will begin to slide. Then, when 
there is relative motion betw^een the drop and the inclined plate, the velocity 
and pressure values inside the drop are numerically computed. The assumptions 
are that the flow within the drop is steady, incompressible and of low Reynolds 
number. The main purpose of studying the patterns of velocity and pressure 
inside the drop is to estimate the coefficient of friction between the flat plate and 
the drop. A proposal for obtaining pressure and velocity values for high Reynolds 
number flows, for the treatment of drop cluster and coalescence, and for obtaining 
unsteady temperature distribution inside the drop is also given. 

1.1 Literature Review 

A survey of the available literature has been presented below as per the following 
sections: (1)- Threshold for movement of a drop down an inclined flat plate, (2) 
Droplet growth and coalescence. 

1.1.1 Threshold for movement of a drop down an inclined 
flat plate 

Furmidge (1961) studied the sliding of liquid drops on inclined solid surfaces 
and also gave a theory for spray retention. He did the study for spray liquids. 
He concluded that surface properties of the spray liquid/solid combination are 
among the most important factors controlling the retention of spray liquids on 
solid surfaces. The spray liquid he used was water and solutions of surfactants, 
while solid surfaces he used was wax and cellulose acetate surfaces. A theory, 
supported by experimental results, has been evolved by him in which he tried 
to explain the movement of drops in terms of the size of the drop, the angle of 
tilt of the surface, the air /liquid surface tension, and the advancing and receding 
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contact angles. Thus for a given liquid on a given surface. 

mgsina _ ^ 
w 


where H is given by 


H = 7az,(cosi9h - cos(9.4) 


( 1 . 1 ) 

( 1 . 2 ) 


For a particular liquid/solid combination H may not be a constant. 
m is the mass of the drop. 

a is the angle of tilt of the surface necessary to produce sliding of the drop. 
w is width of the drop on the solid surface. 

Dimitrakopoulos and Higdon (1999) studied the yield conditions for the grav- 
itational displacement of three-dimensional fluid droplets from inclined solid sur- 
faces through a series of numerical computations. The study considers both sessile 
and pendant droplets and includes interfacial forces with constant surface ten- 
sion. Their study seeks the optimal shape of the contact line which yields the 
maximum displacing force (B'j = B(jsin/5), B^^ is the Bond number and 0 is the 
angle of inclination. This maximum displacing force is that for which a droplet 
can adhere to the surface. 

Kim et al. (2002) have studied the sliding velocity of a liquid drop which 
partially wets a solid surface. This drop will slide down when the plane is tilted 
beyond a critical inclination. Experiments for measuring the steady sliding veloc- 
ity of different liquids of drops were performed by them. Then they constructed 
a scaling law which predicts the sliding velocity given the physical properties, 
wetting characteristics, and size of the drop. When the sliding velocity is low and 
drop distortion due to inclination is small, the scaling law is shown to correctly 
model the functional dependency of the measured sliding velocity. 


1.1.2 Droplet growth and coalescence 

The formation of a distribution of various size droplets is a characteristic feature 
of many systems from thin Aims to fog and clouds. Family et al (1989) have 
presented the results of their investigations of the kinetics of droplet growth and 
coalescence. In general, droplet formation occurs either by spontaneous nucle- 
ation or by growth from heterogeneously distributed nucleation centers, such as 
impurities. Two models have been introduced by them to describe these two 



1.2 Scope of the Present Work 


5 


types of processes. In the homogeneous nucleation model droplets can form and 
grow anywhere in the system. The authors also introduce a heterogeneous nu- 
cleation model for studying processes in which droplets only form and grow at 
certain nucleation centers which are initially chosen at random. Simulations, scal- 
ing theory, and a kinetic equation approach for describing the two models have 
been presented by them. The theoretical predictions are found to be in excellent 
agreement with the simulations. 

When two drops of radius R touch, surface tension drives an initially singular 
motion which joins them into a bigger drop with smaller surface area. The motion 
is always viscously dominated at early times. Eggers et al. (1999) have focused 
on the early-time behavior of the radius of the small bridge between the two 
drops. They have also studied numerically the case of coalescence with an external 
viscous fluid analytically and the case of equal viscosities. 

1.2 Scope of the Present Work 

Earlier work primarily deals with generation of drop shapes above and below 
the inclined and horizontal flat plates. It has not included the stability analysis 
for a drop kept on an inclined plate or hanging below from an inclined plate. 
Modeling of flow inside the drop sliding down the flat inclined plate has also not 
been reported. 

The present work is concerned with not only the generation of drop shapes 
above and below the inclined and horizontal plates, but also involves their stability 
analysis. Moreover, when the drop is in motion, detailed modeling (for steady, 
incompressible and low Reynolds number flow) is carried for obtaining velocity 
and pressure values inside the drop. The coefficient of friction is also estimated 
for the sliding motion of the drop. 

A proposal for studying flow fields at high Reynolds numbers, solution of the 
unsteady energy equation, and treatment of a drop cluster and coalescence is also 
included in the thesis. 
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1.3 Objectives 

The specific objectives of the present study are as follows: 

1. The stability analysis for a sessile or pendant drop, kept on an inclined flat 
surface. 

2. Solution for pressure and velocity in circular and deformed drops for Stokes 
flow. 

3. Unsteady temperature distribution for the deformed drop. 

4. The coefficient of friction for circular and deformed drops. 

1.4 Thesis Organization 

The present thesis has been organized in the following manner: 

1. Chapter 1 presents introduction, literature review and scope of the present 
research. 

2. Chapter 2 describes the determination of drop shapes on an inclined surface 
under static conditions, and stability analysis 

3. Chapter 3 describes mathematical modeling of flow in a sliding semi-circular 
drop. 

4. Chapter 4 presents numerical solution for flow in a sliding semi-circular 
drop. 

5. Chapter 5 describes the grid generation methodology for a drop of a de- 
formed shape. 

6. Chapter 6 presents the derivation of pressure- velocity equations in trans- 
formed coordinates and boundary conditions. 

7. Chapter 7 describes the numerical solution of the governing equations in 
transformed coordinates. 
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8. Chapter 8 is a proposal for studying flow fields at higher Reynolds num- 
bers. solution of the energy equation, treatment of a drop cluster and drop 
coalescence. 

9. Chapter 9 reports conclusions and scope for future work. 


Chapter 2 

Determination of drop shapes on 
an inclined surface under static 
conditions 


Much work has been done on generating drop shapes on inclined surfaces. Ex- 
tensive modeling and simulation has been reported for the relevant physical phe- 
nomena in static and moving drops. 

2.1 Survey of approaches reported in the liter- 
ature 

Much literature is available on the shape of the drops on different configurations 
of the solid surface. Orr et al. (1975) determined interface shapes by using the 
surface representation z — z{x,y) in Cartesian coordinates. Their representa- 
tion was unsatisfactory because the function becomes multi-valued at certain tilt 
angle, specificall when the advancing angle exceeds 90°. Multiple values should 
be avoided because they require the solution of coupled mathematical problems 
of considerably greater complexity. Brown et al. (1980) suggested the spheri- 
cal coordinate system r = f{d,4)) to solve the Young-Laplace equation. They 
implemented a Galerkin finite element method to solve for the shape of a three- 
dimensional drop on an inclined plate. Drop shapes for different volumes and 
different tilt angles of the plate were simulated. It was observed that smaller 
drops deform less than larger ones, keeping other aspects same. Karmakar (2001) 
has also generated drop shapes both above and below a flat inclined surface. He 
generated the drop shapes on the inclined surface and then transformed them to 
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the horizontal. 

2.2 Formulation 

For generating static drop shapes, one has to implement the Laplace- Young equa- 
tion. Let us consider a certain control volume that is occupied entirely by a fluid 
and is fixed in space. As the fluid flows, molecules enter and leave the control 
volume from all sides carrying momentum and thus imparting to the fluid that 
occupies the control volume, at a particular instant in time, a normal force. Fur- 
thermore, short-range intermolecular forces cause the molecules that are located 
on either side of the boundary of the control volume to be attracted, thereby 
generating an effective frictional force. The force exerted on an infinitesimal sur- 
face element of the boundary of the control volume is called the surface stress 
or traction^ and will be denoted by f. Clearly, the value of f will depend upon 
both the position and orientation of the infinitesimal surface element (Pozrikidis 
1997). 

1 



Figure 2.1; Schematic illustration of an interface between two fluids. 

In general, the tractions exerted on the two sides of an interface between two 
fluids labeled 1 and 2 have two different values, with a corresponding discontinuity 


Af = f(i)-f(2) = (aW-(r(2)).n 


( 2 . 1 ) 
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where n is the unit normal vector pointing into fluid 1, as shown in Figure 2.1. and 
and are the stress tensors within the two fluids evaluated at the interface. 
The direction and magnitude of the discontinuity of the interfacial traction Af 
depend upon the mechanical properties of the interface, which are determined 
by the physico-chemical properties of the fluids and molecular structure of the 
interface, and are therefore afiFected by the presence of surface-active substance es. 
An equation that relates Af to the velocity field, the properties of the fluids, and 
the shape and thermodynamic properties of the interface, is called a constitutive 
equation for the discontinuity in the interfacial traction. 

The most common type of interfacial behavior pertains to uncontaminated 
interfaces between two immiscible fluids characterized by isotropic surface tension 
7 , which may be regarded as a kind of energy per unit surface area or surface 
pressure. The physical origin of surface tension may be traced to differences in 
the attraction forces between the molecules of the two liquids. In general, the 
surface tension decreases as the temperature is raised, and vanishes when the 
temperature reaches the boiling point. 

To derive the constitutive equation for Af, we assume that the interfacial 
stratum has negligible mass and write a force balance over a small section of th e 
interface D that is bounded by the contour C, as shown in Figure 2.2, requiring 
(Pozrikidis 1997), 


f AfdS+ f jt 
J D J C 


X ndl = 0 


( 2 . 2 ) 


where n is the unit vector normal to D pointing into 1, t is the unit vector 
tangential to C, and t x n = b is the binormal vector as shown in Figure 2.2. 

Furthermore, the domain of definition of the surface tension and normal vec- 
tor can be extended from the interface into the whole space. The first extension 
can be done in an unrestricted manner, whereas t.he second extension is done by 
setting n = VF/ | VF |, where the equation F{x’^y, 2 ) = 0 describes the instan- 
taneous location of the interface. A variation of Stokes’s theorem states that for 
any arbitrary vector function F, 

f Fxtdl= [ (nV • F - (VF) • n)dS 
Jc Jd 


(2.3) 
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Figure 2.2: Surface tension acting along the edges of a test section of an interface. 


Applying Equation (2.3) with F = 7n and combining the resulting expression 
with Equation (2.2), we obtain 

[ Aids = / [nV • (711) - (V(7n)) • n]d5 (2.4) 

J D J D 

We now take the limit as D shrinks down to a point, expand out the derivatives 
within the integrand on the right-hand side, and thus obtain the desired consti- 
tutive equation 


Af = 7 nV • n — (n X V 7 ) x n 


(2.5) 


The first and second terms in the right hand side of Equation (2.4) express, 
respectively, discontinuities in the normal and tangential directions. When the 
surface tension is uniform, the tangential component vanishes, and the jump in 
the interfacial traction points in the normal direction. 

The divergence of the normal vector on the right-hand side of Equation (2.4) 
is equal to twice the mean curvature of the interface, denoted by 

V • n = 2Km (2.6) 

By definition, the mean curvature is positive when the interface has a spherical 
shape with 2 lying inside and the normal vector pointing outward. 
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The normal component of the traction and thus the pressure undergo a dis- 
continuity across an interface. To compute the jump in pressure in terms of the 
velocity, we resolve Af into its normal and tangential components, and then we 
can have an expression for normal component as 


Af = [—pi -f- P 2 + 2n • (piVu^^^ — paVu*^^) • n]ii -f n x Af x n (2.7) 

Projecting Equation (2.7) onto the normal vector n produces the jump in inter- 
facial pressure in terms of the normal component of Af and the viscous normal 
stresses, 


Ap = Pi — p 2 = — Af • n -f 2n • (pi Vu^’'^ — ) ■ n (2.8) 

where Af is given by an appropriate interfacial constitutive equation. 

For example, when the fluids are stationary or inviscid, the second and third 
terms on the right-hand side of Equation (2.8) are absent. Assuming that the 
interface exhibits uniform isotropic tension, and using Equation (2.5), we find 
Ap = 2,Km'f 

The Navier-Stokes equation in an inertial frame of reference for a fluid that is 
stationary or translates with a uniform velocity u = U(t) takes the simplified form 

cnj „ 

p— = - Vp + pg (2.9) 

Assuming that the density of the fluid is uniform and solving for the pressure, we 
obtain (Pozrikidis 1997) 

P = P (g- • x + A(t) (2.10) 

where A(t) is a time-dependent constant whose value must be found by requiring 
an appropriate boundary condition. 

Using the definition, Equation (2.1), the pressure distribution given in Equa- 
tion (2.10) with vanishing acceleration, and the constitutive equation for the jump 
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in the interfacial traction given in Equation (2.5). we write 

Af = • n = (p 2 ~ Pi)i^ 


( 2 . 11 ) 


where Ap = P 2 — pi- Ai and .42 are two constants, n is the unit normal vec- 
tor pointing into 1, and is the mean curvature of the interface. Rearranging 
Equation (2.11), we obtain the Laplace- Young equation 

Ap 

2«m = — g-x-f-R (2.12) 

7 


where B = (.42 — Ai)/')' is a new constant with dimensions of inverse length. 

For drop on an inclined surface the symmetricity of the drops is destroyed. 
Hence we need a two dimensional formulation that is applicable for unsymmetric 
interfaces. For a two dimensional surface the mean curvature is given by 


1 f" 

2(l + /'2)l 


(2.13) 


Using the Laplace- Young equation (2.12) 


1 f" 

2(l-p//2)| 


7 


g-x + B 


(2.14) 


This on simplification gives 

/"= (1+/'*)^ (2.1o) 

Substituting the capillary number I = {'yJApgf^'^ in Equation (2.15), the gov- 
erning equation is derived as 

r = (I - b) (1 + f ? (2.16) 

To solve Equation (2.16) numerically, we introduce the parameter ‘0, which is 
defined by tan0 = and ranges from zero at the centerline of the drop to 
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the value of the contact angle q at the contact line. Performing the following 
calculations 


and 


f" = 

dx dv dx 
d(— tan^’) dip 
dip dx 


1 dv 
cos^vj dx 


dy dy dx _ dx 

dip dx dip 


(2.17) 


(2.18) 


Substituting Equation (2.17) and Equation (2.18) in Equation (2.16), we get 

I dip (y 


cos^ lb dx 


= (p-B) (1+tanV)’ 


(2.19) 


or 


1 dip 


cos^ip dx 
dip 
dx 
dx 
dip 


\r / cos^yj 
\P J cosip 


cosip 


From Equation (2.18) and Equation (2.20), we get 

dip 

Hence, the first order system is derived as 


dx 

dip 

dip 

where w 


( 

1-B) 

!.20) 

, we get 

sin'0 


-B) 

1 as 



cosip 


w 


sinip 


w 


( 2 . 20 ) 


( 2 . 21 ) 


( 2 . 22 ) 
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2.3 Earlier Method 


Equation (2.22) is to be solved with the boundary conditions 

xiaA) = xo 
y{aA} = 0 

where tp = to — O/e (2.23) 


aA is advancing angle and is receding angle 
The volume constraint is also applicable, namely 

Volume of the drop 


/ 


x'^dy = 


TT 


x[aA) = xq 


(2.24) 


where xq is the starting x coordinate of the drop as when kept on the horizon- 
tal plate. To solve the Equation (2.22), the following procedure was used by 
Karmakar (2001). 

1. There are several unknowns. One is the radius of the drop and second is 
the constant B. In addition and are also unknown. To simplify the 
problem we assume that the contact line is pinned and does not change 
with inclination. Then we solve for zero inclination to get the radius of the 
drop. 

2. Start the integration with the initial condition, Equation (2.23), and an 
assumed value of B. Integrate the equation using the order Runge- 
Kutta method till ip = a. Then we calculate the volume. If it matches the 
prescribed volume, terminate the iterations; otherwise change the value of 
B using the Newton-Raphson method. 

3. After getting the value of B and radius, we start the integration with the 
initial condition ip = —a. Continue the integration till the angle equiv- 
alent to the plate inclination is reached. Reject all the points ahead and 
continue the integration till the curve (representing the drop) intersects the 
inclination line. Stop the integration. 

4. Transform the inclined drop in normal plane. To satisfy the conditions. 
Equation (2.24), we need scaling. To satisfy the radius constraint we per- 
formed X scaling and to conserve the volume, the y coordinates were scaled. 
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2.4 Method Employed in the Present Work 

Equation (2.22) is to be solved with the boundary conditions 

a : ( 0 . 4 ) = xq 

y{aA) = 0 

where tp = Oi a to an (2.25) 

a A is advancing angle and an is receding angle 


Oi-A = 

an — Op Oq 


(2.26) 


The volume constraint is also applicable, namely 



Yblume of the drop 

TT 


i(a^) = lo 


(2.27) 


where xq is the starting x coordinate of the drop as when kept on the horizontal 
plate. The numerical algorithm is as follows 


1. There are several unknowns. These are namely, constant B, a a, Op, xq, 
volume. These are provided as starting parameters for the problem, except 
constant B. To simplify the problem we assume that the contact line is 
pinned and does not change wuth inclination. 

2. Start the integration with the initial condition, Equation (2.25), and an 
assumed value of B. Integrate the equation using the 4*^ order Runge- 
Kutta method for xp = aA to ij) — an- Then we calculate the volume. If it 
matches the prescribed volume, terminate the iterations; otherwise change 
the value of B using the Newton-Raphson method. 

3. After getting the value of B, we start the integration with the initial con- 
dition xp = aA to xp = an- Stop the integration. 
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Karmakar (2001) generated the drop shapes above and below the inclined 
surface and transformed them to the horizontal. On the other hand the method 
employed in the present work, generates the drop shapes on the incline for differ- 
ent angles of inclination of the surface. This method gives a more clear picture of 
the drop shape on or below an inclined surface than method followed by Karmakar 
(2001). Also in experimental work the results can be easily related to the drop 
on or below the incline rather with the drop shape which has been transformed 
to the horizontal. 

2.5 Numerical Technique 

The numerical technique employed in solving the first order coupled ordinary dif- 
ferential equations, Equation (2.22) is the 4^^ order Runge-Kutta method. For a 
single dependent variable, we have 


dt 

yitn) 


Vn 


(2.28) 


The 4^^ order Runge-Kutta algorithm for solving Equation (2.28) is 


!/n+l = v- + >‘[j + j + j + j)+0(h'‘) 
h = f{tn,yn) 

^2 = /(^n + pil/n + ^■^) 

h = f{in + -^yyn + h-^) 
h = /(in + h,yn + hkz) 


(2.29) 


For Equation (2.22), there are two dependent variables x and y, and only one 
independent variable tp. 
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2.6 Formulation for drop shape below a surface 

The Laplace- Young equation for a pendant drop is 

1 f" 




(2.30) 


when the drop is below an inclined flat plate 


/" 

3 - = — -9y + B 


2(l + f2)| 


(2.31) 


This on simplification gives 

/" = (-ypy-s) (1 + /'")" (2.32) 

Substituting the capillary number I = ( 7 /Apy)^^^ in Equation (2.27), the gov- 
erning equation is derived as 

r=(-|-B)(l + /'')^ (2.33) 

To solve Equation (2.28) numerically, we introduce the parameter tp, which is 
defined by tan^* = — and ranges from zero at the centerline of the drop to 
the value of the contact angle a at the contact line. Performing the following 
calculations 


fif _ ^ 

dx d'lp dx 

_ d{—ta.nip) dip _ 1 dip 

dip dx cos^ip dx 


(2.34) 


and 


dy _dy dx _ dx 

dip ~ dx # “ 


(2.35) 
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Substituting Equation (2.29) and Equation (2.30) in Equation (2.28). we get 


1 dv 
cos^v dx 


(1 + tanV)^ 


(2.36) 


or 


1 dp 
cos'^p dx 
dp 
dx 
dx 
dp 



From Equation (2.30) and Equation (2.32), we get 

dy _ sinp 

dp 

Hence, the first order system is derived as 


(2.37) 


(2.38) 


dx 

cosp 

dp 

w 

dy _ 

sinp 

dp 

w 

where w = 



2.7 Earlier Method 

Equation (2.39) is to be solved with the boundary conditions 


x{aA) = xo 
y{aA) = 0 

where p = q-a to — (2.40) 
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0,4 is advancing angle and aR is receding angle The volume constraint is also 
applicable, namely 


J x'^dy 

x{aA) 


Volume of the drop 

TT 

Xo 


(2.41) 


where xq is the starting x coordinate of the drop as when kept on the horizon- 
tal plate. To solve the Equation (2.39), the following procedure was used by 
Karmakar (2001). 

1. There are several unknowns. One is the radius of the drop and second is 
the constant B. In addition Q !.4 and cur are also unknown. To simplify the 
problem we assume that the contact line is pinned and does not change 
with inclination. Then we solve for zero inclination to get the radius of the 
drop. 

2. Start the integration with the initial condition. Equation (2.40), and an 
assumed value of B. Integrate the equation using the 4*'* order Runge- 
Kutta method till ip = a. Then we calculate the volume. If it matches the 
prescribed volume, terminate the iterations; otherwise change the value of 
B using the Newton-Raphson method. 

3. After getting the value of B and radius, we start the integration with the 
initial condition ip = -a. Continue the integration till the angle equiv- 
alent to the plate inclination is reached. Reject all the points ahead and 
continue the integration till the curve (representing the drop) intersects the 
inclination line. Stop the integration. 

4. Transform the inclined drop in normal plane. To satisfy the conditions 
(2.41) we need scaling. To satisfy the radius constraint we performed x 
scaling and to conserve the volume, the y coordinates were scaled. 
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2.8 Method Employed in the Present Work 

Equation (2.39) is to be solved with the boundary conditions 

x(aR) = Xo 
y{aR) = 0 

where = auto a a (2.42) 


aA is advancing angle and an is receding angle 


a A = 9p — 9 c 

OiR = 9c 


(2.43) 


The volume constraint is also applicable, namely 

/ „2 j.. _ Volume of the drop 
X dy — 

TT 

x{aR) = Xo (2.44) 

where xq is the starting x coordinate of the drop as when kept on the horizontal 
plate. The numerical algorithm is as follows. 

1. There are several unknowns. These are namely, constant B, aA, 9c, 9p, Xo, 
volume. These are provided as starting parameters for the problem, except 
constant B. To simplify the problem we assume that the contact line is 
pinned and does not change with inclination. 

2. Start the integration with the initial condition, Equations (2.40), and an 
assumed value of B. Integrate the equation using the 4*^ order Runge- 
Kutta method for 'ip = aR to ip = aA- Then we calculate the volume. If it 
matches the prescribed volume, terminate the iterations; otherwise change 
the value of B using the Newton-Raphson method. 

3. After getting the value of B, we start the integration with the initial con- 
dition Ip =■ an to Ip = aA- Stop the integration. 
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2.9 Numerical Technique 

Numerical technique adopted for solving Equation (2.39) is 4^^ order Runge- 
Kutta method of Section 2.5. 

2.10 Stability Analysis 

Properties of fluids used are as follows: 


Tablel. Properties of Fluids used. 


S.No. 

Fluid used 

Density (kg/m^) 

Surface Tension (N/m) 

Capillary Number (m) 

1 

Water 

1000 



2 

Bismuth 

10401.8 

0.39557 

0.0019699 


Mercury 

13545.9 

0.486226 

0.0019138 


Stability analysis for the sessile drop above the incline and pendant drop below 
the incline, are separately carried out. 

2.10.1 Sessile Drop 

Center of gravity is calculated by taking the mass of any element as pdxdy. Then 
center of gravity coordinates {xcg,ycg) are given by 


Xcg — 
Vcg — 


f f xpdxdy 
f J pdxdy 
J fypdxdy 
J f pdxdy 


(2.45) 


The advancing angle is calculated from the slope of the line joining the first and 
the second coordinate of the equation describing the drop shape. The receding 
angle calculation is done from the slope of the line joining the last and the last 
but one coordinate of the equation describing the drop shape. 

Moment calculations for forces acting on a sessile drop, as shown in Figure 
2.3, are given below: 
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Normal Reactio^T 



Figure 2.3; Force analysis for a sessile drop. 


Normal Reaction = mgcosdp + 2BRi + 2^fRi[sin{9^4,£)) + sin(6RE)] (2.46) 


where mg is the magnitude of the weight force of the liquid, comprising the drop, 
B is the excess pressure term (difference between the ambient and pressure inside 
the drop), Ri is the radius of the span which the drop covers on the flat inclined 
plate, 7 is the surface tension, 6 ad is the advancing angle of the drop, Ore is the 
receding angle of the drop. 

Taking moment about the advancing angle of the drop, we have 


Moment = (Normal Reaction)i?i + mgsin6p[ycgCos6p— \ Xi — x^g \ sindp] 

-mg(xi - Xcg) - BRi - ijR^sinBRE (2.47) 


where xi is the x coordinate of the advancing angle point of the drop. 

On the other hand, considering only the forces which make the drop move along 
the incline, moment is calculated as follows: 


Momentl = mgsin$p[ycgCos6p- 1 Xi - Xcg 1 sin$p] - A'yRlsinORE (2-48) 

For a three dimensional drop, only moment of those forces which make the 
drop move down the incline, as shown in Figure 2.4, is calculated as 
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Figure 2.4: Base of a three dimensional drop representing forces, which make the 
drop move down the incline. 


Moment2 


mgsindp[ycgCosdp— \ Xi — Xcg\ s 
—47^2^ / (1 — cos(j))sin Q. 

J<t>=o 


Up] 


0AD + 9 RE . 

h * 


dd) (2.49) 


where (j) is the base angle (assuming to be a circle), ranges from 0 to tt. Also h is 
the total no of steps in which 6ad+Gre is divided and i goes from 1 to (/i+ 1)/2, 
if h is odd 'whereas i goes from 1 to {h/2) + 1 or 1 to (h/2) — 1, when h is even. 


2.10.2 Pendant Drop 

For pendant drop case (drop below a flat inclined surface), as shown in Figure 
2.5, the stability calculations are as follows. The center of gravity is calculated 
by taking the mass of any element as pdxdy. Then center of gravity coordinates 
{Xcg, Vcg) are given by 
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_ f f^pdxdy 
1 1 Pdxdy 
f f y pdxdy 
1 1 pdx:dy 


(2.50) 


The advancing angle is calculated from the slope of the line joining the last 
and the last but one coordinate of the equation describing the drop shape. The 
receding angle calculation is done from the slope of the line joining the first and 
the second coordinate of the equation describing the drop shape. 



Figure 2.5: Force analysis for a pendant drop. 


Moment calculations are given below: 


Normal Reaction = 2BRi — mgcosBp + 2'yRi[sin{9AD) + sin{diiE)] (2.51) 


where mg is the magnitude of the weight force of the liquid comprising the drop, 
B is the excess pressure term, Ri is the radius of the span which the drop covers 
on the flat inclined plate, 7 is the surface tension force, 9ad is the advancing 
angle of the drop, 9re is the receding angle of the drop. 

Taking moment about the advancing angle of the drop, we have 


= (Normal Reaction)ili + mgsin9p[ycgCos9p- j Xi — Xcg 1 sin9p] 
+mg I xi - Xcg I -BRi - 4'yRlsin9RE (2.52) 


Moment 
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On the other hand, considering only the forces which make the drop move along 
the incline, moment is calculated as follows: 


Moment 1 


mgsin6p[ycgCosdp— | Xi — Xcg \ sindp] + mg I Xi — Xcg | —BRi 
—A^fF^sinOnE (2.53) 


For a three dimensional drop, only moment of those forces which make the drop 
move down the incline, as shown in Figure 2.3, is calculated as 


Moment2 


mgsindp[ycgCos6p— \ xi — Xcg \ sinOp] + mg | Xj — Xcg 
— 47i2f / (1 — cos<j))sin 

J (2>=0 


AD 


d.AD + SrE . 

-X 


h 


dd) 


-BRi 

(2.54) 


where <j) is the base angle (assuming to be a circle), ranges from 0 to tt. Also h is 
the total no of steps in which Oad + Ore is divided and i goes from 1 to (h + 1)/2, 
if h is odd whereas i goes from 1 to (/i/2) + 1 or 1 to {h/2) — 1, when h is even. 

2.11 Results and Discussion 

Figures 2.6 and 2.7 show sessile and pendant drops respectively, generated by 
using axysymmetric formulation in the Laplace- Young equation. These drop 
shapes are for comparison with the two dimensional unsymmetric formulation for 
0° plate inclination, shown in Figures 2.8 and 2.9. As in both two dimensional 
formulation and axysymmetric formulation, for 0° plate inclination, the drop 
shape should be the same. But as the plate inclination is different from 0°, two 
dimensional formulation should be used to generate the drop shapes. Figure 2.6 
is for 3 different liquids namely, mercury, bismuth and water, having contact 
angles 120°, 135° and 60° respectively, with 2 different volumes 25 and 50mm°. 
Similarly, Figure 2.7 is for the same parameters, except that bismuth and mercury 
volumes are 4 and 5mm°, as for Figure 2.6, but former is for sessile drop case and 
latter is for pendant drop case. 
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Figure 2.8 shows the change in sessile drop shapes for 3 different liquids 
namely, mercury, bismuth and water, having volumes 25mm^. 50mm^. 50mm^ 
and contact angles 120°. 135°. 60° respectively. It shows the change in drop 
shapes as plate inclination is changed from 0 to 30°. As already discussed, the 
advancing angle and receding angle are Qc. — and 6c with respect to the inclined 
plate respectively. So, as the plate inclination increases, the advancing angle 
should decrease and receding angle should remain constant, with respect to the 
inclined plate. 

Figure 2.9 shows the change in pendant drop shapes for 3 different liquids 
namely, mercury, bismuth and water, having volumes 25mm°, 50mm°, 50mm° 
and contact angles 120°, 135°, 60° respectively. It shows the change in drop 
shapes as plate inclination is changed from 0 to 30°. As already discussed, the 
advancing angle and receding angle are dc and dc — dp with respect to the inclined 
plate respectively. So, as the plate inclination increases, the advancing angle 
should remain constant and receding angle should decrease, with respect to the 
inclined plate. 

Figures 2.10, 2.11, 2.12 shows for mercury, bismuth, water sessile drops respec- 
tively, variation of advancing angle, receding angle, moment and momentl with 
respect to the plate inclination. Generally, advancing angle increases, receding 
angle decreases, moment and momentl increases with respect to plate inclination. 
Here the drop shapes have been simulated with 6c and dc — 6p as advancing angle 
and receding angle with respect to horizontal, but w'hen calculated with respect 
to the incline flat plate they are 6c — 6p and 6c respectively. 

Figures 2.13, 2.14, 2.15 shows for Mercury, Bismuth, Water pendant drops 
respectively, variation of advancing angle, receding angle, moment and momentl 
with respect to the plate inclination. Generally, advancing angle increases, re- 
ceding angle decreases, moment and momentl increases with respect to plate 
inclination. Here the drop shapes have been simulated with dc — dp and dc as ad- 
vancing angle and receding angle with respect to horizontal, but when calculated 
with respect to the incline flat plate they are dc and dc — dp respectively. 



D(op Height (m) 
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(c) 




(d) 



Figure 2.6: For zero plate inclination and axysymmetric formulation (a) mercury- 
drop of volume 25mm^ and contact angle 120° (b) mercury drop of volume 50mm° 
and contact angle 120° (c) bismuth drop of volume 25mm° and contact angle 
135° (d) bismuth drop of volume 50mm° and contact angle 135° (e) water drop 
of volume 25mm° and contact angle 60° (f) water drop of volume 50mm° and 
contact angle 60°. 
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(a) (b) 







Figure 2.7: For zero plate inclination and axysymmetric formulation (a) mercury 
drop of volume 4mm® and contact angle 120°, (b) mercury drop of volume 5mm® 
and contact angle 120°, (c) bismuth drop of volume 4mm® and contact angle 
135°, (d) bismuth drop of volume 5ram® and contact angle 135°, (e) water drop 
of volume 25mm® and contact angle 60°, (f) water drop of volume 50mm® and 
contact angle 60°. 
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Figure 2.8: Sessile drop shapes for diflferent plate inclinations from 0-30° of (a) 
mercury of volume 25mm^ and contact angle 120°, (b) bismuth of volume 50mm^ 
and contact angle 135°, (c) water of volume 50mm° and contact angle 60°. 




2.11 Results and Discussion 


31 




-5E-05 " V 


0 “ 

10 “ 

20 “ 

30“ 




0.0002 0.0004 


Figure 2.9: Pendant drop shapes for different plate inclinations from 0-30° of (a) 
mercury of volume 25mm° and contact angle 120°, (b) bismuth of volume 50mm° 
and contact angle 135°, (c) water of volume 50mm° and contact angle 60°. 
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Figure 2.10: mercury sessile drop of volume 25mm^, contact angle 120° (a) varia- 
tion of advancing angle with respect to plate inclination, (b) variation of receding 
angle with respect to plate inclination, (c) variation of moment with respect to 
plate inclination, (d) variation of momentl with respect to plate inclination. 
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Figure 2.11: bismuth sessile drop of volume 50mm^, contact angle 135° (a) varia- 
tion of advancing angle with respect to plate inclination, (b) variation of receding 
angle with respect to plate inclination, (c) variation of moment with respect to 
plate inclination, (d) variation of momentl with respect to plate inclination. 
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Figure 2.13: mercury pendant drop of volume 25mm®, contact angle 120° (a) vari- 
ation of advancing angle with respect to plate inclination, (b) variation of receding 
angle with respect to plate inclination, (c) variation of moment with respect to 
plate inclination, (d) variation of momentl with respect to plate inclination. 
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Figure 2.14: bismuth pendant drop of volume 50mm^, contact angle 135° (a) vari- 
ation of advancing angle with respect to plate inclination, (b) variation of receding 
angle with respect to plate inclination, (c) variation of moment with respect to 
plate inclination, (d) variation of momentl with respect to plate inclination. 
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Figure 2.15: water pendant drop of volume 50mm^, contact angle 60° (a) variation 
of advancing angle with respect to plate inclination, (b) variation of receding angle 
with respect to plate inclination, (c) variation of moment with respect to plate 
inclination, (d) variation of momentl with respect to plate inclination. 



Chapter 3 

Mathematical modeling of flow in 
a sliding semi-circular drop 


The governing equations, namely continuity, momentum, and pressure equations 
for a semi-circular drop along with boundary conditions are presented below 
(Muralidhar et at 1995). 

3.1 Equations of Fluid Motion in a semi-circular 
region 

The governing equations for fluid flow in r-9 coordinates are given below. The 
mass balance equation is 


dp d{(M*) pu* I d{pv*) 

dP dr'*' r* r* dO 


(3.1) 


where u*, and v* are the components of velocity in r* and 6 direction respectively. 
The superscript ★ indicates that the respective quantities are in dimensional form. 
The momentum equations in component form are 
r* Component 


P 






dp* 

dr* 


+ p 


\7^u* 






(3.2) 


where X*, Y* are the components of body force in r* and 6 direction respectively, 
p is viscosity and p* is pressure. 
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9 Component 


dv* 

IF 


u*v* 


= Y* 


1 dp* 
r* dO 






du* 

~d6 


(3.3) 


where 


du* _ du* ^ ^ / du* \ V* ( du* \ 

dt* dt* ^ y Qj-* j r* \ d6 ) 

2 * d^u* 1 f du*\ 1 f d‘^u*\ 

V ^ ^ J^F^ \ W ) 


(3.4) 

(3..5) 


The following assumptions have been additionally enforced in the governing 
equations. The flow is of low Reynolds number (Re< 1), steady, and incompress- 
ible. The body forces namely X* and 1'* are neglected. The validity of these 
assumptions is based on the fact that the drop is physically small in size. For a 
small drop, surface forces are dominant in comparison to body forces. The mass 
balance equation, momentum equation, and pressure equation along with these 
assumptions and the boundary conditions for semi-circular drop placed on a flat 
horizontal surface are given below. 

The mass balance equation in r* 9 coordinates is 


du* u* 1 / dv* \ _ Q 

dr* r* r* \ 50 J 


(3.6) 


The momentum equations in component form are 
r* Component 



(3.8) 
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From the mass balance and momentum equations, we obtain the governing equa- 
tion for pressure as 


d'^p* 1 / dp'’ \ 1 / d'^p* \ _ Q 

dr*^ J-* J 7-*2 ^ Q02 J 


(3.9) 


3.2 Dimensionless Governing Equations in the 
Physical Domain 


Dimensionless governing equations in r and 9 co-ordinates are given below, where 
superscript * denotes the dimensional quantities, the non-dimensional compo- 
nents of velocity are u = u*/uo, v = v*/uo where uq is the plate velocity. 

The non-dimensional radius, pressure are r = r*/ri, p = p* — pl^i,/ pu^, where 
ri denotes the maximum radius of a data point in the physical domain, p is the 
density of fluid of the drop. Re = pu^ri/p, where Re is Reynolds number and p 
is the viscosity of the liquid phase of the drop. 

The mass balance equation in dimensionless form is 


du u I 

« — ' ^ — 

or r r 



= 0 


The momentum equations in component form are given below^ 
r Component 


dp 1 d^u 1 1 f d'^u\ 2 ( 

dr Re dr"^ r ySr/ 



(3.10) 


(3.11) 


9 Component 

1 f dp\ 1 I ^ dv\ 1 /'df^v'\ 2 du\ 

^ r / "*” Re dr^ r \9r/ \d9^ ) \d9 J 

The boundary conditions for velocity are 
For ^ = 0 


V 


(3.12) 


u = l,v = 0 


(3.13) 


For 9 — IT 


u = — l,u = 0 


(3.14) 



3.2 Dimensionless Governing Equations in the Physical Domain 41 


At fluid interface 


u-n = 0 


(3.15) 


and 


5(u • t) 
dn 


0 


(3.16) 


where u is velocity vector and n is outward drawn unit normal vector on the drop 
surface and t is the unit tangent vector on the drop surface. Let F(x. y) = 0 be 
the equation describing the drop shape, then n =yF/ 

The pressure governing equation is 


^4. If 4.L( = 0 

dr ‘2 r \dr J \d9'^ J 


(3.17) 


Pressure boundary conditions are: 
At the fluid interface 


p = 2km''/ 


(3.18) 


where 7 is surface tension and km is mean curvature of the drop in physical 
domain. Let F{x,y) = 0 ov y = f{x) be the equation describing the drop shape. 
Then km is defined as (Pozrikidis 1997) 


km — 


f" 


2(1 + 


(3.19) 


If F{r, 0) = 0 be the equation of the drop, then km for a semi-circular geometry 
is 


km — 


1 J_ 

2ri 


(3.20) 


For 6 = 0 and 0 = tt (at the plate), the pressure boundary conditions are 


dp _ 2r 'd‘^v 1 f 1 / ^ f 

d9 Re dr"^ r \ 5r / ”^ \ d9‘^ ) \d9 J 


(3.21) 


As V is zero and u is constant at the plate 9 = 0 and 0 = tt, the derivatives of v 


and u with respect to r are zero. 



Chapter 4 

Numerical solution of the Flow 
Equations for semi-circular drop 


The dimensionless form of pressure governing equation is 


^ r \dr J I 66^ _ 


(4.1) 


The boundary condition for pressure are 
For 0 = 0 and 0 = tt (at the plate): 


dp _2r \d'^v 1/^^ 1 fd‘^v\ 2 fdu\ 

^ ^ Re r Vdrj r2 


For r = 0 


p = average of all values at neighboring points 


For r = 1 ( from Equation (3.18) ) 



(4.2) 


(4.3) 


(4.4) 


where We is Weber number and is expressed as puln/^. 
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For solving the governing equation for pressure along with boundary condi- 
tions, we adopt the finite difference technique. The radial distance of the cylin- 
drical coordinates is discretized in to n — 1 steps, with the index of discretization 
in radial direction being i. The angular length of tt radians of the cylindrical 
co-ordinates is discretized in to m — 1 steps, with the index of discretization in 
angular direction being j. There are m grid points in the radial direction and n 
grid points in the angular direction. 

Equation (4.1) can be written in symbolic form as 


4 ^ 

5r2 


dp ^ d^p 


= 0 


(4.5) 


Discretizing this canonical form Ave get 


‘^Phj Pi-lJ 

+ B 


Pi,j+1 ~ ^PiJ 4" PiJ-1 

2Ar 

= 0 

{Adf 



(4.6) 


or 


pihj) 


2[A{Mf -b C'(Ar)^] 

(A0)^ 




(2.4 4- SAr)pi+i j+ 


Pi-1, - BLr) + C'(Ar) (pij+i +pij_i) 


(4.7) 


where 


A = 
B = 


1 

1 



C 


(4.8) 
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Discretizing the boundarj" conditions, we get 
For ^ = 0 


Pi,l = Pi.2 - A 


V aF 



B[Uj,2 — 


(4.9) 


For 6 = T[ 


Pi,m — Pi,m-1 + -4 


^i,m 






A9 




m-lj 


(4.10) 


For r = 0 


p = average of all values at neighboring points (4-11) 


For r = 1 


1 


We 


(4.12) 


Dimensionless form of the r component of the momentum equation is 


1 d'^u 1 f 9u\ 1 ^ d^u\ 2 / 9v\ u' 

Re r \9r / '*' \d9‘^ J \d9 ) r^_ 


dp 

dr 


(4.13) 


The boundary conditions for r component of velocity are 
For 0 = 0 


u = 1 


(4.14) 


For 0 = TT 


u = — 1 


(4.15) 
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For r = 0 


u — linear interpolation between -1 and 1 


(4.16) 


For r = 1 


u • n = 0 


(4.17) 


Equation (4.13) can be written as 


.d'^u -du ^d'^u 





(4.18) 


where 


A 

B 

C 

D 

E 

F 


1 

1 




(4.19) 



4- /? 


1 

< 

1 

“T -LJ 

2Ar 


Discretizing Equation (4.18) by finite differences we get 
A 

+C 

L 

tt> [ Pi+1 J ~ Pi-1 J 

Wr 


^ij+1 

4- n 

1 

rH 

1 

1 

r-l 

+ 

i 

1 

> 

to 

i 

-T JL/ 

Ae j 


+ Euij 


(4.20) 
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or 

1 

~ [2A{Aef + 2C{Arf - E(Aef{Arf] 

Ui+ij^^(2.4 + BAr) + (ui.i,j^^^(2A - BAr)4- 

A0 

-D--^(Ar) {Vij^i — 

F^(AS)"fe«j-pi.ij) (4.21) 

Discretized boundary conditions are 
For 0 = 0 

Ui,x = 1 (4.22) 

For 0 = TT 

'^i,m ~ 1 (4.23) 



The dimensionless form of the 0 component of momentum equation, which is the 
0 component of velocity governing equation, is 


1 1 ^d^v\ 2 _ 1 

Re dr"^ "** r \5r J \50^/ / 


(4.26) 
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The boundary conditions for 6 component of velocity are 
For ^ = 0 


r = 0 


(4.27) 


For 9 = Tr 


v = 0 


(4.28) 


For r = 0 


u = 0 


(4.29) 


For r = 1 


5(u • t) 

dn 


= 0 


(4.30) 


Equation (4.26) can be written as 


where 


.d^v ^dv d‘^v du _ _ 


dp 

m 


A 

B 

C 

D 

E 


1 

1 




Til 


hJ 


fjf*. .2 

1 

. .2 
‘hJ 

Re 


’t,3 


(4.31) 


F 


(4.32) 
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Discretizing Equation (4.31) by finite differences, we get 


.4 

C 




-hF 

i-'i-fij j 

. 

(Ar)^ 


2Ar j 




+ D 


- 

(AS)' 

_ 

2A0 

F 

PiJ-l 

2A0 





- Evij 


(4.33) 


or 




1 

{2A{Adf + 2C'(Ar)^ + E(Ar)^(A^)^) 

- Di\e)+ 

«'iJ+l^(2C + DA«l) - F^iArfipij^, - ftj-i) 


(4.34) 


The discretized boundary conditions are 
For 0 = 0 


Vi,i = 0 


(4.35) 


For 0 = TT 




= 0 


(4.36) 


For r = 0 


vij = 0 


(4.37) 
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For r = 1 


^n.j — 




(4.38) 


4.1 Iteration Algorithm 

1. Grid points are generated in the physical domain. The pressure and mo- 
mentum equations are coupled and are to be iteratively solved. 

2. Initial guesses for u, v and p are taken to be zero. 

3. Pressure equation is solved for zero initial values of u and v. 

4. Numerical algorithm used for solving p, u and v is the Gauss-Seidel iterative 
method. Convergence criteria used is such that the ratio of the difference 
of the present and previous value, to the average value at each grid point is 
less than 10“®. 

5. Then u is solved with the present values of p and v. 

6. Numerical algorithm and convergence criteria used for u is as in Step 4. 

7. Next, V is solved with the present values of p and u. 

8. Numerical algorithm and convergence criteria used for v is as in Step 4. 

9. But as the governing equation for pressure contains the velocity terms, the 
velocity components should possess a unique value. Therefore the procedure 
from 5-8 is repeated until the convergence criteria for u and v at every grid 
point is met. 

10. Next the coefficient of friction is calculated. It is defined as the ratio of 
average frictional force to the average normal force. 

11. After we get a fixed value of u, v for a particular value of p at each grid 
point, these values of velocity components now serve as inputs for the next 
iteration for solving pressure. 
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12. Steps 5-10 are again repeated. 

13. The above algorithm continues until convergence criterion for coefficient of 
friction is met. Convergence criteria for co-efficient of friction is as gi\'en in 
Step 4. 

14. Lastly, we get the converged values of p, u, v and the coefficient of friction. 

4.2 Code Validation 


The code validation for Laplace equation in a semi-circular region has been carried 
out. The model problem has Dirichlet boundary conditions. 

The governing differential equation is taken as 


d'^p ^^dp d^p _ 


(4.39) 


where 


A = 1 



C = A (4-«) 


with the boundary conditions are 

For ^ = 0 

Pi,i = 0 

(4.41) 

For 6 = TT 

Pi,m ~ 0 

(4.42) 

For r = 0 

Phj = 0 

(4.43) 

For r = 1 

Pn,j — 1 

(4.44) 

The analytical solution of the model problem is 

Pij - - f ^sm(tij) -h sm(3tij) -f- ^ sin(5tij)^ 

(4.45) 


The comparison of the numerical and analytical solutions is shown in Figure 
4.1. The analytical and numerical results are similar. The dotted line shows the 
analytical solution and the solid line shows the numerical solution. 



4.3 Grid Independence 


o 



Angle (radians) 


Figure 4.1: Code validation of Numerical Solution against the Analytical Solution. 

4.3 Grid Independence 

For 9 = 72°, the profiles for pressure for all radii are compared for grids of 
Meaningful results 101 x 101 and 51 x 51. The results can be seen on the 51 x 51 
grid in Figure 4.2. The nature of pressure values for 101x101 and 51x51 are 
same, hence all the calculations have been done for 51x51. 
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e 51 X51 



Radius 


Figure 4.2: Grid independence test for 51x51 and 101x101. 


4.4 Coefficient of Friction 

The coefficient of friction expressed in terms of dimensional quantities is 

^ ‘i’-* + ^ {^) '*>■* 


coefficient of friction = 


P(f^l)5+(P*-Pamb)2?'l 


(4.46) 


When expressed in terms of non-dimensional form we get 



4.4 Coefficient of Friction 


53 


coefficient of friction = 


MjZl^li^ldriuori) 4- u Ll-i dr (i) driupr-j 


P(§r-f)g~2km-(2ri) 


i:4.4i 


where 


2fc„ = V-(^) (4.48) 

where F(x, y) = 0 is the equation of curve in Cartesian coordinates. If F(r, 0) = 0 
is the equation of the semi-circle then, 

2km = - (4.49) 

ri 


Therefore, the expression for coefficient of friction for a semi-circular drop is 


coefficient of friction = 


^Ir=o d: (i) ^ (ID driuon) 

p(frf)5+ i7(2ri) 

(4.50) 


or 


where 


/ imr^g \ , ( 27puori \ 

2jutio ) ~ V ) 


(4.51) 


+ /r=-l H w) 

^ iw to) 


(4.52) 


Key 


Pr = 


We = 


- 

/i 

Up 

\/(pn) 

puln 


7 


(4.53) 
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For lOmm^ drop and Wq =lm/s, Fr is 35.7434. 

The table for coefficient of friction for a semi-circular drop, with We and Re var- 
ied over three orders of magnitudes, can be seen on the next page. The coefficient 
of friction increases with decreasing Reynolds number because of an increase in 
viscosity. It decreases with decreasing Weber number (increasing surface tension) 
because of an increase in the bulk pressure of the drop. 


Table2. Variation of Coefficient of Friction with Re and We. 


S.No. 

We 

Re 

Coefficient of Friction ! 

1 

0.01 

0.01 

2.789 i 

2 

0.1 

0.01 

27.961 

3 

1.0 

0.01 

279.457 

4 

0.01 

0.1 

0.291 

5 

0.1 

0.1 

2.789 

6 

1.0 

0.1 

27.946 

7 

0.01 

1.0 

0.003 

8 

0.1 

1.0 

0.291 

9 

1.0 

1.0 

2.787 


4.5 Results and Discussion 

Figure 4.3 shows pressure contours and velocity vectors in a semi-circular drop. 
The boundary has a constant pressure. Inside the drop, the pressure must be pos- 
itive, except where it is negative. This negative pattern are due to recirculation. 
The boundary condition for velocity is that there is no flow in the direction of 
the normal vector to the bounadry and there is no shear stress at the boundary. 
The velocity vectors are computed with respect to the plate. Similarly Figures 
4.2 to 4.9 have the same boundary condition for pressure and velocity, except in 
all these figures Weber number and Reynolds number vary over three orders of 
magnitude. 
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Figure 4.5: For We=1.0 and Re=0.01 (a) pressure contours (b) velocity vectors 
relative to the plate. 


(a) (b) 




Figure 4.6: For We=0.01 and Re=0.1 (a) pressure contours (b) velocity vectors 
relative to the plate. 
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(a) (b) 




Figure 4.7: For We=0.1 and Re=0.1 (a) pressure contours (b) velocity vectors 
relative to the plate. 


(a) 



(b) 


1.5 


Figure 4.8: For We=1.0 and Re=0.1 (a) pressure contours (b) velocity vectors 
relative to the plate. 
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(a) 




Figure 4.9: For We=0.01 and Re=1.0 (a) pressure contours (b) velocity vectors 
relative to the plate. 



Figure 4.10: For We=0.1 and Re=1.0 (a) pressure contours (b) velocity vectors 
relative to the plate. 


Y Coordinate 



Figure 4.11: For We=1.0 and Re=1.0 (a) pressure contours (b) velocity vectors 
relative to the plate. 



Chapter 5 

Grid Generation 


The grid generation procedure that transforms the deformed drop geometry to 
a semi-circular geometry is described in the present chapter. Figure 5.1 shows 
mapping of physical domain (deformed drop) to computational domain (semi- 
circular region). Let the body-fitting coordinate transformation be of the form 


i 



Figure 5.1: Physical domain (deformed drop) is mapped to computational domain 
(semi-circular region) 

^ = ^(r, 0) and rj = r]{r, 9). From chain rule 



^rdr + ^gd9 
7]rdr + r]ed9 


& 6 1 Mr 1 

Vr Ve \ \ d9 ) 


(5.1) 


where (^r 5 ^ 0)5 iVri'He) are partial derivatives of f and 77 with respect to r and 9 
respectively. In a similar way, considering r = rj) and 9 = 9{^, 77 ), we obtain 



1 / 

9^ e,j \ \ dr] } 


(5.2) 


Combining Equations (5.1) and (5.2) 
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{- 1 = 


- ^n]\ & fg 1 


e. ^ 9,1 _ 


nr 


1 / 1 
ns \\ d9 j 


'■5.3 j 


or 


Cr 6 

Vt ne 


rr, 

9^ 9, 


J_ [ 1 

L j 


(-5.4) 


where | J| = r^9, — 9^r,, is the determinant of the Jacobian matrix 


r p — 

[Ji - 


9^ 9, J 


(5.5) 


For one-to-one mapping to exist, the determinant I Jj should be finite and non- 
zero. From Equation (5.4), it is clear that the partial derivatives of ^ and rj are 
related to those of r and 9 through the relations 


?r — 9,/ I J| ; — T,/ I Jj : 

r^r = -9j\J\: n9 = n/\J\. 


(5.6) 


5.1 Grid Generation Technique 

For regular geometries, the grid generation can be done anahdically through suit- 
able algebraic transformation. For irregular geometries, obtaining transformation 
functions is difficult; also, there is no guarantee that the particular transforma- 
tion employed may produce an acceptable grid. A more general procedure for the 
algebraic generation of grids is based on interpolation. Here, the coordinates of 
the interior nodes are obtained by interpolation between the prescribed boundary 
data. A very useful technique for this purpose is the Transfinite Interpolation. 
As given by Muralidhar et al. (1995), the general algorithm of transfinite inter- 
polation can be stated as; 



5.1 Grid Generation Technique 


62 


1. Place grid points on the boundaries of the domain as desired. 

2. Apply unidirectional interpolation in f-direction for ^-direction) betv.-een 
the boundary grid data given on the curves ^ = 0 and ^ = 1 (or r? = 0 and 
77 = 1) and obtain the coordinates for every interior and boundary point. 

3. Calculate the mismatch between the interpolated and the actual coordinates 
on the 77 = 0 and 77 = 1 (or ^ = 0 and C = 1) boundaries. 

4. Linearly interpolate the difference in the boundary point coordinates and 
find the correction to be applied to the coordinates of every interior point. 

5. For each point, final coordinate value is equal to uni-directionally interpo- 
lated value - correction. 

The transfinite interpolation is a very simple and powerful numerical procedure. 
The main weakness of this procedure is that the slope discontinuities at the 
boundary propagate to the interior and spoil grid smoothness. However, since the 
metric derivatives are also evaluated numerically, the effect of a slope discontinuity 
is not severe, as in the case of analytical transformation. 

Figure 5.2 shows the mapping procedure for grid points from physical domain 
(deformed drop) to computational domain (semi-circular region). 



Figure 5.2: Mapping procedure for grid points from physical domain (deformed 
drop) to computational domain (semi-circular region) 

This procedure is described as follows: 


T] = 9 


(5.7) 


Then ^ = ^{r, 6 ) is determined for each 9 . 

If P is any grid point in the physical domain and P' is the grid point where line 
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joining origin (0) to P when extended meets the boundary. Q is the correspond- 
ing grid point on the line joining origin ( O) to Q. OQ is given by 


0Q = 


OP 

op- 


(5.8; 


Chapter 6 

Governing equations in 
transformed coordinates 


Dimensionless form of the governing pressure equation in transformed coordinates 
can be derived to be of the form: 


ae drj^ d^drj 




(6.1) 


where 


A 

B 

E 

C 

D1 


Sr ' ^9 

2 f^r»7r ■ 


r2 ) 


r 

„ , r)ee 

Tjrr 1 77 


( 6 . 2 ) 


There boundary conditions in ^-t] coordinates are 
For 7] = 0 and 77 = tt 
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where 


.4 

B 

C 

D 

E 

F 

G 


el 


4 


,^eriO 


iee 

J.2 

m 

r2 

2^ 

r2 

.Ve 


(6.4) 


For ^ = 0 

p = Average of pressure values at surrounding grid points (6. .5) 


For ^ = 1 


P = 


We 


( 6 . 6 ) 


Dimensionless form of the governing equation for r component of velocity is 


^ d^u du du , r^^v di' 

5^ OT] 


(6.7) 
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where 


-4 

B 

C 

D 

E 

F 

G 

H 

AA 

BB 


Cl 


Sr 


2 (^^rVr 


,c 

sr 


^ee 


) 


T r‘ 

. Vr . me 



Re(Cr) 

Re(77r) 


The boundary conditions are 
For ry = 0 


Ui,i = 1 


For rj = TT 


'^i,m — 1 


(6.8) 


(6.9) 


( 6 . 10 ) 


For e = 0 


Uij = U2,l + (j - 1) 


U2,m - ^ 2,1 


m — 1 


( 6 . 11 ) 
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For ^ = 1 


u 


n.j 



■ n ,;-: - \ 

2A77 ) 


{ 6 . 12,1 


The dimensionless form of the governing equation for 6 component of velocity is 


d‘^v d'^v d‘^v 
^ ar?2 ■ dr]d^ 


, dv dv du du 


Hv 


(6.13) 


where 


.4 

B 

C 

D 

E 

F 

G 

H 

Col 

Col 



y.1 


2 



+ - -r ^ 

„ , '7r , r]6Q 

?7rr d i 5“ 

r 





(6.14) 


The boundary conditions are 
For 77 = 0 


Vi,i = 0 


(6.15) 
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For ^ = 0 


For C = 1 



Simplifying all equations for the case where r) = 6 and only fy {or r^) is non-zero, 
we obtain the following equations: 

Dimensionless form of the governing pressure equation in transformed coordinates 
can be derived to be of the form: 


de didri 



(6.19) 
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where 


.4 = 0 



C = ^ 

T 

D1 = 0 (6.20) 


There boundary' conditions in ^-rj coordinates are 
For 77 = 0 and 77 = tt 

^ = il ( A~ R— ■ r-^^D—^E—\ - f n 
dr] Re \ 5 ^^ dvjd^ dr]) ^ 


( 6 . 21 ) 


where 


.4 = 0 



D = 0 
£■ = 0 
F = 0 



(6.22) 


For ^ = 0 

p = Average of pressure values at surrounding grid points (6.23) 


For e = 1 


P = 


1 

We 


(6.24) 
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Dimensionless form of the governing equation for r component of velocity is 


d^u d^u 


■<■<1 


BB 


^ d^u du du dv . dv 
^ d^drj ' dr] ' drj 

dr. 


(6:2o) 


where 


.4 

B 

C 

D 

E 

F 

G 

H 

.4.4 

BB 


= e 

j_ 

r2 

= 0 

— §L 
r 

= 0 
= 0 

J.2 

“ 

= Re(^r) 
= 0 


(6.26) 


The boundary conditions are 
For 77 = 0 


Ui,l = 1 


(6.27) 


For 77 = TT 


tii,m — 1 


(6.28) 
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For ^ = 0 


u 


■j 


= U2.1 ^ (j - 1 J 


U-7 


■TC - ^ 2.1 

m — 1 


For ^ = 1 


‘■n.j 


_ 


rn.j +1 - rn.j-l 

2Ae 


i6.29) 


(6.30) 


The dimensionless form of the governing equation for d component of velocity is 




(6.31) 


where 



(6.32) 


The boundary conditions are 
For ?7 = 0 


For 77 = 


For C = 0 


For ^ = 1 


'^i.m — 0 


Vij = 0 


n/ [ 1 + 


Un-lJ-lTn-l.j ~rn-l,j -2 

2 rn-.i 


2rn-i,j^e 

V + V 2r„_i,jAe ) 


+ ^-’n-lj-l 


v/ i + (TsSrri?^y 


t^n,7^nj+l ^n.j— I 

2r„jA0 



Chapter 7 

Numerical solution of the 
governing equations in 
transformed coordinates 


Dimensionless form of the governing pressure equation in transformed coordi- 
nates, with i being the index of variation in the radial direction (r or ^). both in 
physical and computational domains, and j being the index of variation in the 
angular direction {9 or rj), both in physical and computational domains. 


4 ^ 

9^2 


d^p d^p 

d^dp 




(7.1) 


where 


E = + 


c 

D1 



(7.2) 
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scretizing the governing differential equation for pressure, we get 
.4 


Pi-ij - \Pi.j-i - ‘2pi.j ^Pi.j-C\ 


V — J 


iAlJ} 


Pi-rl,j^l Pi-rlJ-l Pi-l,j-l ~V Pi-l,j-l | ^ | Pi-'-j Pi-'^-j i 


-D1 


Pi,j~rl ~ Pi.j-- '\ _ « 

2A?7 j 


2AC 


(T.3J 


PiJ 


1 

[2.4(A?7)^ + 2S(Ae)^] 

+ CA?) + + D1A^)+ 

- -D1A^)+ 

(A0(A/7) , „ 1 

^ iPt+lj-t-l -Px+lj-l -Pi-lj+l +Pi-lj-i; 


Pi-ij— ^(2.4 - C Af)+ 


(7.4) 


e boundary conditions are 
r 77 = 0 
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Pi , 2 
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Tje [Re 


dp 
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^ + E-\ 

dC dn) 


(7.5) 


ere .4, B, C, B, B, F, G have been defined in Equation (6.4). 

' Tj — ir 
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Re 
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dp 


(7.6) 




For ^ = 0 


p = Average of pressure values at surrounding grid points 


For ^ = 1 


^ We 


Dimensionless form of the governing equation for r component of velocity is 


where 
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(7.10) 


76 


Discretizing by finite difference, we get 


or 
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The boundary conditions are 
For rj = 0 


Ui,l = 1 


(7.13) 
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For /] = rr 


^j.m — 1 


(7.14; 


For ^ = 0 


Uij — U 2,1 -r {j — I 


^2.m — « 2.1 
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For ^ = 1 
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Dimensionless form of the governing equation for d component of velocity is 
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Discretizing by finite differences, we get 
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The boundary conditions are 
For 77 = 0 


t‘j,i — 0 


(7.21) 


For T) = 7t 


^'i,m — 0 


(7.22) 


For ^ = 0 


vij = 0 


(7.23) 


For ^ = 1 
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7.1 Coefficient of Friction 


Coefficient of friction for the deformed drop can be calculated as 


coefficient of friction = 





(7.25) 
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where 


Rev = 


Fr = 


We 


pupri 

P 

Up 

s/im) 

1 


(7.26) 


For lOmm^ drop and no=lm/s. Fr is 35.7434. 

Below is table for coefEcient of friction for a deformed drop, with We and Re 
varied over three orders of magnitudes. The coefficient of friction increases with 
decreasing Reynolds number because of an increase in viscosity. It decreases with 
decreasing Weber number (increasing surface tension) because of an increase in 
the bulk pressure of the drop. 

Tables. Variation of Coefficient of Friction with Re and We. 


S.No. 


Re 

Coefficient of Friction 

1 

0.01 

0.01 

1.912 

2 

0.1 

0.01 

19.441 ^ 

3 

1.0 

0.01 

201.124 

4 

0.01 

0.1 

0.272 

5 

0.1 

0.1 

1.931 

6 

1.0 

0.1 

19.434 

7 

0.01 

1.0 

0.019 

8 

0.1 

1.0 

0.288 

9 

1.0 

1.0 

1.924 


It can be noted that drop deformation does not have a major influence on coeffi- 
cient of friction, when compared to that for a semi-circular drop. 


7.2 Results and Discussion 

Figure 7.1 shows pressure contours and velocity vectors in semi-circular region 
(computational domain) and deformed drop (physical domain), for Weber number 
0.01 and Reynolds number 0.01. The boundary has a constant pressure. Inside 
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the drop, the pressure must be greater than that on the boundary. The boundary 
condition for velocity is that there is no flow in the direction of the normal vector 
to the boundary and there is no shear stress at the boundary. The velocity vectors 
are computed with respect to the plate. Similarly Figures 7.2 to 7.9 show pressure 
contours and velocity vectors in semi-circular region (computational domain) and 
deformed drop (physical domain), for the same boundary condition for pressure 
and velocity, except in all these figures Weber number and Reynolds number \’ary 
over three orders of magnitude. As Weber number increases from 0.01 to 1.0. the 
values of pressure contours decreases. With increasing Reynolds number, velocity 
is increased. 
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(a) 



(c) 



(b) 



(d) 

2 - 



Figure 7.1: Pressure contours and velocity vectors (relative to the plate) for 
We=0.01 and Re=0.01. (a), (c): computational domain; (b), (d): physical do- 
main. 
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Figure 7.4: Pressure contours and velocity vectors (relative to the plate) for 
We=0.01 and Re=0.1. (a), (c): computational domain; (b), (d): physical domain. 
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the plate) for 
vsical domain. 
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Figure 7.8: Pressure contours and velocity vectors (relative to the plate) for 
We=0.1 and Re=1.0. (a), (c): computational domain; (b), (d): physical domain. 
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Chapter 8 

Proposal for studying flow flelds 
at high Reynolds numbers and 
unsteady Temperature 
distribution 


The stream function-vorticity formulation is a useful approach for solving two- 
dimensional incompressible flows. One of the main attractive features of this 
method is that it does not involve the solution of the pressure field. In the 
stream function-vorticity method, pressure is completely eliminated by cross- 
differentiation of the momentum equations. 

8.1 Solution for flow flelds at high Reynolds num- 
ber 

The mass balance equation is 


u du 1 dv 
r dr r 50 


( 8 . 1 ) 


r momentum equation is 


du vdu _ dp 1 d^u 1^ _ — 

dr r do r dr Re dr^ r dr r^ dO^ r^ dO 


( 8 . 2 ) 
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6 momentum equation is 


dv , udv uv _ Idp 1 'd^v 1 dr , 1 d^r 2 r] , 

^ dr r dd r r 89 He dr^ rdr~^ r'^dO'^ ' r^J 


To obtain the stream function- vorticity {p — uj} equation, p is eliminated from 
Equation (8.2) and Equation (8.3) by differentiating Equation (8.2) with respect 
to 9 and Equation (8.3) with respect to r and subtracting one from the other. 
The resulting equation is expressed with vorticity ui as the dependent variable 
which is defined by 


d'^ip Idif 1 d^p 
dj -2 r dr r'^ 89^ 


(8.4) 


The result is 


du vduj _ 1 'd^o; Ido; 1 d^<j 

^ dr r 89 Re _ dr^ r dr ^ 89“^ 


where the radial and tangential velocity components are 


Idp 
r 89 


and V = — 


dp 

dr 


The boundary conditions for the velocity are 
For 0 = 0 


(8.5) 


( 8 . 6 ) 


u = l and r = 0 


(8.7) 


For 9 = IT 


u = —1 and r = 0 


( 8 . 8 ) 
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For r = 0 


u = linear interpolation between -1 and 1 and r = 0 (8.9) 


At the interface 


u - n = 0 and 


d(u • t) 
dn 


= 0 


(S.IO) 


From these boundary conditions for velocity components, the boundary condition 
for the stream function can be found. 

Transforming the Equation (8.4) into computational domain 


4^ 

de 


drf' d^rj drj 


= 0 


where 


^ = e4 

^ + ^ 

C = 2 

D = + | + ^ 

E = Vrr + ^ + '^ 


( 8 . 11 ) 


( 8 . 12 ) 


Transforming the Equation (8.5) into computational domain 
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where 
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C 

D 

E 




Vr 


S77 


2 &Vr + 


^rr 

TjfT 


r2 

^ee 


+ 


r 

nr 


- Re (u^r - 

‘Hee f , V \ 
-^-Re(uvr^-r]e) 


(8.14) 


8.2 Solution of Pressure Equation 

In the stream function-vorticity method, to obtain pressure at each grid point 
for viscous flow, it is necessar}’’ to solve an additional equation for pressure. This 
equation is derived by differentiating Equation (8.2), r momentum equation, with 
respect to r and Equation (8.3), d momentum equation, with respect to 9. Com- 
bining with Equation (8.1), mass balance equation, we finally get the governing 
equation for pressure. The boundary conditions can also be found in terms of 
stream function. Thus pressure can be solved for each grid point, knowing the 
value of stream function. 

8.3 Solution of Unsteady Thermal Energy Equa- 
tion 

Dimensionless form of the unsteady thermal energy equation is 

dT dT IdT 1 fd'^T IdT 1 d^T\ 

It de ~ ?e\dr^ ^ r dr ^ 89^ ) ^ 


where u, T, t are the dimensionless form of velocity, temperature difference and 
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time scaled by Uo, AT and ri/uo respectively. 

Initial condition is 
At t=0 

r = l 18.16) 


Boundary conditions are 
For ^ = 0 and ^ = 77 

r = 0 (8.17) 


For r = 0 


r = o 


(8.18) 


For r = 1 


r = 1 


(8.19) 


Transformation to the computational domain leads to the following dimensionless 
form of the governing equation in computational domain 


dt 




where 


A 

B 

C 

D 
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Pe 

1 

Pe 
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Pe 

1 

Pe 

1 

Pe 




C ^ 


Srr -r _ ^ 


{nr 
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) - - y 

r T- J 
nr . n96\ 

— + -V - 'U77r - u 
r / 
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r 
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r 


( 8 . 20 ) 


( 8 . 21 ) 
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Discretizing Equation (8.16) by finite differences, we get 

Ti/ = 1 (8.22) 


Discretizing the boundary conditions by finite differences, we get 
For ?7 = 0 


Ti/ = 0 


(8.23) 


where i goes from 2 to n — 1 
For r} = -K 


T- P = n 


(8.24) 


where i goes from 2 to n — 1 
For ^ = 0 


Ti/ = 0 


(8.25) 


where j goes from 1 to m 
For ^ = 1 


T’ .p _ 1 


(8.26) 


where j goes from 1 to m, p denotes particular time step. Discretization variable 
in radial direction (^) is i, its value ranges from 1 to n. Discretization variable in 
angular direction ( 77 ) is j, its value ranges from 1 to m. 

There are three popular schemes available for solving initial value unsteady prob- 
lems. These are 

1. Euler method (or explicit method) 


2. Crank-Nicholson method 
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3. Pure implicit method 

Of these three methods, Crank-Nicholson is the most accurate and pure implicit 
is the most stable. To make an optimum choice out of these three, pure implicit 
method is usually preferred over the other two. In the pure implicit method, the 
time derivative at the new' time is used to move ahead in time. Thus 

BT 

r"'‘ = r»+— I (8.27) 

Here is substituted from Equation (8.20) in Equation {8.27). hence we 

get an expression in which can then be solved iteratively using the Gauss- 
Seidel method. 

8.4 Drop Cluster and Coalescence 

When a single drop is surrounded by many drops, then there is much probability 
that this drop will coalesce with one or more of its surrounding drops. Thus in 
the study of multiple drops, the phenomena of coalescence becomes significant. 

When two drops of radius R touch surface tension drives an initially singular 
motion which joins them into a bigger drop with smaller surface area. Eggers et 
al. (1999) showed that this motion is always viscously dominated at early times. 
They focused on the early-time behavior of the radius of the small bridge 
between the two drops. The flow is driven by a highly curved meniscus of length 
27rr„i and width A around the bridge, from w'hich it can be concluded that 
the leading-order problem is asymptotically equivalent to its tw'O-dimensional 
counterpart. For the case of inviscid surroundings, an exact two-dimensional 
solution (Hopper, 1990) shows that A oc and r,„ ~ {t'y/Trr])ln[tj/{r]R)]\ and 
the same is true in three dimensions. (Hopper, 1990) also studied the case of 
coalescence with an external viscous fluid. A significantly different structure is 
found in which the outer-fluid forms a toroidal bubble of radius A oc at 
the meniscus and ~ it'y/4Trr])ln[t'y/{rjR)]. This basic difference is due to the 
presence of the outer-fluid viscosity, however small. 

When a droplet of radius n touches or overlaps a droplet of radius r 2 , a new 
droplet is formed, centered on the center of mass of the tw'O original drops. If this 
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droplet overlaps one or more other droplets, they are coalesced and this procedure 
continues until no overlaps remain. 


8.5 Results and Discussion 

Results for the unsteady temperature problem are shown in Figure 8.1. Figure 
(a), (b), (c), (d) shows the temperature contours (Pe=0.01) in a deformed drop at 
11^^, 23^^, 35^^, 47^^ time steps and Figure (e) shows temperature contours at 
47^^ time step in the computational domain. Figure (f) shows the temperature 
evolution for a particular grid point. Steady state was reached in 47^^ time step. 



Figure 8.1: Temperature contours (Pe=0.01) in a deformed drop (a) at 11^^ time 
step, (b) at 23^^ time step, (c) at 35^^ time step, (d) at 47^^ time step, (e) in the 
computational domain (semi-circular region) at 47*^^ time step, (f) temperature 
evolution for a particular grid point with time. 




Chapter 9 


Conclusions and Scope for Future 
Work 


Drop shapes above and below an inclined and horizontal flat plate, are generated. 
Stability analysis of sessile and pendant drops were carried out for different plate 
inclinations. It was observed that with increasing plate inclination, the drop 
becomes more unstable. Three liquids were considered namely, mercury, bismuth 
and water, besides having different densities and coefficient of surface tension, 
their contact angles are also different. The plate inclination was varied from 0 to 
30°. A formulation for the moment generated by forces that slide the drop down 
the incline is also presented. 

Numerical solution for pressure, velocity are obtained for a semi-circular drop. 
The results are reported for three orders of magnitude variation, of Weber number 
and Reynolds number. For a deformed drop, the grids are mapped from the phys- 
ical domain to computational domain. The governing equations and boundary 
conditions are then transformed in the computational domain. After obtaining 
the solution for pressure and velocity in the computational domain, then using 
transformation relations, solution for pressure and velocity in the physical do- 
main is obtained. Unsteady temperature contours have also been obtained for a 
deformed drop. 

The numerical solutions obtained are for a two dimensional geometry and for 
Stokes flow. A proposal for solving pressure and velocity for Reynolds number 
greater than one, and drop coalescence has also been provided. 

Results show that the coefficient of friction increases as Reynolds number and 
Weber number decreases. The effect of deformation is seen in the local pressure 



100 


and velocity fields, but the overall impact on the coefficient of friction is small. 

Scope for future work consists of obtaining pressure and velocity fields for 
Reynolds number greater than unity. Pressure and velocity solution for three 
dimensional drop is also recommended for future work. Modeling for coalescence 
is required both for two dimensional and three dimensional drops. 
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